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Решена задача о 2)-распределении биополимеров при изоэлектрическом фокусировании в прямоугольной 
электрофоретической камере в иммобилизованном градиенте рН. Программное обеспечение, созданное на 
языке ЕгееРет-—-, позволяет в наглядной форме исследовать картину изоэлектрического фокусирования во 
времени, а также моделировать и оптимизировать реальный эксперимент. 
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Введение. Изоэлектрическое фокусирование (ИЭФ) является одним из наиболее эффективных 
электрохимических методов, применяемых для идентификации компонентов и фракционирования 
многокомпонентных биохимических смесей. Этот метод широко используется в биохимии, различ- 
ных отраслях медицины, генетике, фармакологии, пищевой промышленности. 

Основной задачей математического моделирования ИЭФ является создание моделей, по- 
зволяющих повысить разрешающую способность метода, оптимизировать эксперимент в целях 
экономии материальных и временных затрат. Для этого необходимы наглядные модели, позво- 
ляющие исследовать процесс ИЭФ во времени, а также установить взаимосвязь между динамикой 
процесса и параметрами электрофоретической камеры (ЭК). 

На решение перечисленных задач направлены наши исследования. Полученная модель 
основана на классических математических моделях [1, 2] и содержит новый подход к их решению 
в интегрированной среде разработки РгееРет+- [3]. 

Язык программирования ЕгееРет+-+ выбран как оптимальный для решения указанной за- 
дачи по ряду причин. Во-первых, ЕгееРет-+-+ изначально предназначен для численного решения 
дифференциальных уравнений в частных производных методом конечных элементов в 20- 
пространственном случае. Во-вторых, по сравнению с коммерческими пакетами, предназначен- 
ными для решения уравнений в частных производных, ЕгееРет-+- обладает уникальной возмож- 
ностью доступа ко всем внутренним данным и возможностью создания собственных алгоритмов. 
В-третьих, в настоящее время существует развитая теория решения задач математической физи- 
ки средствами языка РгееРет-+-, в частности, разработан обширный математический аппарат для 
решения задач электрофореза как метода разделения многокомпонентных смесей на отдельные 
компоненты при помощи внешнего электрического поля [4]. В-четвертых, язык ЕгееРет-+-+ отно- 
сительно прост в использовании, поскольку запись алгоритмов на нем близка к С++, а эффектив- 
ность кода по скорости близка к оптимальной. Наконец, в-пятых, язык РЕгееРет--+ обладает раз- 
витым интерфейсом, имеет собственный «визуализатор», поддерживающий рисование сетки, 
изолиний конечно-элементарных линий и векторных полей. 

Применение языка РЕгееГет++ в рассматриваемом случае позволило создать программное 
обеспечение, моделирующее поведение сложных систем биополимеров в заданных условиях ИЭФ. 
Физическая постановка задачи. В электрофоретическую камеру, имеющую форму прямо- 
угольника длиной / и высотой й, помещена смесь № биополимеров в виде пятна заданной конфи- 
гурации. В ЭК задан фиксированный градиент рн. 
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Для каждого из биополимеров известны его константы диссоциации рК“, рК®, а также 
коэффициенты миграции п,. Известно также количество биополимеров Д,, помещенных в ЭК. 


Предполагается, что температура внутри ЭК постоянна и равна Т. 
Требуется исследовать динамику концентрации каждого из № биополимеров во времени 
под действием электрического поля: рассчитать неизвестную концентрацию биополимеров 


Е =Ё, (х,у,1) ‚ Е=Ь2,..., М и визуализировать результаты расчетов с помощью линий уровня. 


Математическая постановка задачи. Предполагается, что диссоциация биополимера ПОД НО- 
мером Кописывается уравнениями: 


МН:ВСООН <> МН.ВСООН-Н*”, —МН,ВСООН <> МН,ВСОО`+Н*. 

Молярные концентрации МН:ВСООН, МН,ВСОО`, МН.ВСООН - положительного, от- 
рицательного и «нейтрального» ионов биополимера — обозначим &, &\, &. Аналитическая 
концентрация биополимера определяется формулой: &, = &1 + +&\,. В равновесном состоянии 
концентрации ионов связаны с аналитической концентрацией следующими соотношениями: 

Е = а , = > = , ы = (- о. > ОЙ , 
Н? й КИК" 
КОКОКО?" КОК К®Н+Н?' 








К 
ея 


где а и а^ — степень диссоциации биополимера; Н — концентрация ионов водорода. 


В условиях иммобилизованного (фиксированного) градиента рН (т.е. функции 
Н=Н(х, у) ‚ заданной в области ДР = {х=[0;1], у=[0;1]} и неизменной во времени) концентра- 
ции биополимеров описываются краевой задачей: 

















Ре -а(р У +, (и —0*,)Уф)=0, =1,2,...М (1) 
с. Аф+Уф.Ус=0, (2) 
о= > Их (1 +0 Е, АР (3) 
[2 бы + ри, п(а т ол =0, (4) 

дп —и 

0ф 0ф 

| =0, || =0, 5 
[5 [$ 
Фо =Ф1, Ф,,=$>, (6) 
[ТЕ (у) =М,, (7) 


ур) 
где о — проводимость среды в ЭК; ф -— потенциал электрического поля в ЭК; У = [>= — гра- 
диент функции; ОР, — коэффициенты диффузии, связанные с коэффициентами миграции 
известным уравнением: 
О; =ЕнН,, 
Где == АТ/Р. 


Уравнение (1) есть уравнение массопереноса, полученное на основании уравнения потока 
биополимера и основного уравнения теории переноса. Уравнения (2), (3) выражают закон Ома, 
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записанный в предположении, что электрический ток в растворе создается миграцией биополи- 
меров, а также ионов водорода. Краевое условие (4) отражает факт, что граница ЭК непроницае- 
ма для электролита и, следовательно, поток каждого из биополимеров через границу равен нулю. 
Краевые условия (5) соответствуют предположению о том, что верхняя и нижняя границы облас- 
ти изолированы, а значит, электрический ток через них равен нулю. Условия (6) означают, что 
для ЭК определена разность потенциала между ее левой и правой границами. Наконец, инте- 
гральное условие (7) отражает тот факт, что масса каждого биополимера в течение всего экспе- 
римента остается неизменной и равной М, . 


Следует учесть, что при интегрировании представленной интегродифференциальной 
краевой задачи требуется также использовать распределение в начальный момент времени 
Ее (х,у,0) в некоторой заданной ограниченной области. 


Система уравнений (1)-(7) представляет собой серьезную проблему для численного реше- 
ния, в первую очередь, из-за наличия интегрального условия (7). Решение ее классическими ко- 
нечно-разностными методами требует предварительных преобразований системы, составления 
сложного алгоритма, значительного объема вычислений. 

В то же время, применение интегрированной среды разработки РгееРет--+ к решению за- 
дачи позволяет минимизировать как предварительные преобразования, так и непосредственную 
вычислительную работу. 

Численная реализация задачи на ЕгееРет++. Алгоритм решения рассматриваемой задачи 
на языке РЕгееРет-+-+ включал в себя шесть основных этапов. Рассмотрим кратко каждый из них. 

Параметризация границ области интегрирования. В соответствии со стандартами 
РгееРет-+-+, граница области была обозначена через Г (замкнутый прямоугольный контур, обхо- 
димый против часовой стрелки); при этом нижнюю и верхнюю границы области обозначали 


В (у=0) и 2, (у=1); правую и левую соответственно [. (х =/) и [,(х=0) (рис.1). 


(0, 1) [3 (1, п) 
УР Ра 
(0,0) Г (1,0) 


Рис.1. Область интегрирования р = {(х,у) | [0,/] х [0,#]} 


Аппроксимация временных производных. В соответствии с требованиями языка 
РгееРет-+-+ задача была приведена к слабой (вариационной) формулировке. Однако непосредст- 
венное (без составления алгоритма) решение возможно только для линейных стационарных за- 
дач. Поскольку рассматриваемая задача моделирования ИЭФ достаточно сложна, потребовалось 
пошаговое решение, позволяющее свести исходную задачу к набору линейных краевых задач. 


Поиск решения задачи (1)-(7) осуществляли на интервале времени [0,7 ]. Был задан на- 


бор значений #, =тт, где т — шаг по времени, а также введены обозначения: 


в (ху) = (%,У,1„), „ =тх. 
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С учетом аппроксимации производной по времени конечной разностью задача (1)-(7) бы- 
ла преобразована к виду: 


т-+1 ый 
































= -ам(р, УБЕ (ой 02 Уф") =0, К=1,2,..,М (8) 
(н”= ) К®кК® 
0: = ' о = г 2 ' (9) 
К®К® +К®Н”ч В Ее ) К®К® +К®Н"Н + (н”" ) 
с"Аф" +Уф"Уо" =0, (10) 
№ 
"= Ура (0 +44 Е" +рнН", (11) 
= 
НЫ тик К т == 
р, —— +Ы, А п(а ри о )Уф =0 7 (12) 
дп ыы 
ди ны ди ме 
ф" | = ф, Г Я, — ф> Г (14) 
ег" (ъу)аиау = М, . (15) 
р 


Таким образом, если &, (х,у)=&, (х,у,) и ф"(х,у)=ф(х,у,!,) известны, то для опре- 


т-+1 


деления функции &,” (х,у)=&,(х,уи„.) имеет место стационарная краевая задача (8)-(15). 
Решая ее для т = 0,1,...,и, можно получить приближенное решение нестационарной задачи. Для 
решения задачи используется неявная аппроксимация, т.е. все члены уравнений, не содержащие 
производной по времени, выбираются в точке #,., = (т-+1)т, а аппроксимация производной по 
времени осуществляется в точке 1#, =тт. 


Переход к слабой (вариационной) формулировке задачи. В основе метода конеч- 
ных элементов лежит приведение задачи к слабой (вариационной) формулировке. 

Первоначально был осуществлен переход к слабой формулировке уравнения (10). Для 
этого использовали исходное уравнение, приводящее к уравнению (10): 


("Уф") =0. 
Умножение исходного уравнения на пробную (тестовую) функцию 9, интегрирование по области 
Р споследующим интегрированием по частям привело к уравнению: 


[[с"Уф”Убчкау + фбо" О ЕО. 
ры г дп 
В соответствии с краевыми условиями (13) 


| О ео: 





УВ " 
а значит, 
[| ° р а | бо" 99-0. (16) 
5 ах дд во, бт 
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Из краевых условий (14) следует, что 





у и =ф,, Ч, =ф.. (17) 


В языке РгееРет-+- имеется большой набор быстрых алгоритмов для решения задач, приведен- 
ных к слабой формулировке. В данном случае для решения задачи (16), (17) был выбран стан- 
дартный код с привлечением Ш-алгоритма: 

рго_ет рЬРИИр/!Ы, у, воет = ГИ) = 


ш2А(ТА)( яюта * (ах(рй1) * Ах(») + 4у(р#й1) * 4у(»))) 
+оп(Ё4, ри = рйП) + оп(Е2, ри = р#2). 

Аналогично была получена слабая формулировка уравнений (8). Умножение каждого из 
уравнений (8) на пробную (тестовую) функцию 0, и интегрирование по частям в области Р при- 
вело к уравнению: 

о, Чу (РЕ НН» (ат — Уф") = 


= = ов (Ру +," (ой — о Уф" ухау + 


НЫ 
д 
в Р- МЕЧИ, — = | 


В силу краевого условия (12) а интеграл в правой части последнего уравне- 
ния равен нулю, а значит, уравнение принимает вид: 


[[° 9 


В соответствии с а языка полученные уравнения необходимо просуммировать 
по всем А, в результате чего получим слабую формулировку задачи, состоящей из уравнений (8), 


(12): 
ео [Аа а, 
ит р. 9х дх ду д 


ы Гы (а - 75 : & 2 | у [, 
Для решения задачи (18) также был выбран код с использованием И\-алгоритма (для простоты 
рассмотрен случай № =2): 
ргоМет рЬ4с(с1, с2, 7, у2, ти = т, воет = ГО) 
= шё2А(ТА)( с1* УЕ ерзЙопт1* (@х(с1) * ах(1) + 4у(сТ) * 4у(1))) 
+тЕ2а(Ти)( ти1* | * с1* (ах * ахС7) + а’ф® * 4у(”))) 
-шЕ2а(ТР)( со * 1) 
+те2а(СТИ)( с2* у2 /Ё- ер5Иоп? * (Ах(с2)* 4х(»2)- 4у(с2)* 4у(у2)))\ 
+щтЕ2а(Ти)( ти2 * | *с2* (Ах(рт) * 4х2) + ду) * 4у(»2))) 
-шё2а(ТИ)( со 42/4 * у2 ); 


Стандартные действия. Для параметризации границ области был использован следую- 
щий код: 





к 


ВР руб + ей - оодуох" |в 80: 








(18) 




















х1 = 0; у1= 0; х2 = и]; у2 = 0; х3 = и]; у3 = #1; х4 = 0; у4 = #1; 
Богаег Г1(1 = 0,1) { х = (1-Е) * х1+1* х2; у= (1-0 * У1-1* у2;}; 
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Богаег Г.2(1= 0,1){ х= (1-#*х2+1*х3З; у= (1-* у2-+1* у3; $}; 
Богаег ЁЗ(Е= 0,1) { х= (1-Е) * хЗ+1*х4; у = (1-1) * уз+#* у4; }; 
Богаег Г.4(1 = 0,1) { х = (1-1) * х4+1*хГу= (1-8) * у4+1*у1; }; 

Триангуляция области (т.е. замена ее конечным набором треугольников) осуществлялась 


посредством оператора 
тезй Ти = БиИатезй(ГЛ(10 * п) Г2(5*п)- [3(10*п) + [4(5*п)) 


и и к генерированию сетки (рис.2). 


ими ИУ 
к АЛЛАААДЕКИИ А 2 








 ОЗАУУДУУАУАУ УДАВУ ВЧ УЬЧЬЧ, ых 
У ДУДУАУ КУ БЗРАУДУДУ ЧБ АЗИИ 


Рис.2. Сетка интегрирования 





В качестве конечных элементов были выбраны кусочно-квадратичные непрерывные ко- 
нечные элементы Уй(Тй,Р2) , обеспечивающие одновременно высокую точность и скорость реше- 


НИЯ. 

Организация цикла движения по времени. Особого внимания при создании про- 
граммного обеспечения потребовало создание цикла движения по времени. Цикл был организо- 
ван с помощью стандартного оператора /юх’(т= 0; т < 1000; т- +), на каждом шаге которого 


требовалось обеспечивать условие нормировки. 

При решении задачи необходимо было иметь в виду, что по физическому смыслу концен- 
трация биополимера есть величина неотрицательная. Одним из приемов, который можно исполь- 
зовать для сохранения неотрицательности, является обработка на каждом временном шаге дис- 
кретного набора концентраций в узлах сетки. Если шаг расчета достаточно мал и при вычислении 
решения в каких-либо узлах появляются отрицательные значения концентраций, то их можно 
принять равными нулю. 

Пусть на некотором временном шаге 1, существует решение, которое обозначим 


АЯ У). Определим отклонение среднего значения вычисленной концентрации от действи- 
тельного значения по формуле: 





55 & (жуи,)-М, Ми4у =Е,, 


где М, - масса А-го биополимера; о — мера области О. 


пед 


Концентрацию на (т +1) -м шаге вычисляем по формуле: 


С Е = СЕ (ху, )- =} ” 
Очевидно, что в этом случае масса и а автоматически: 
ДЕ х, у)ахау = К=1,...,М, 


что соответствует условию нормировки (15). 
Перерасчет концентрации осуществлялся с помощью следующих кодов: 
Лог (ий ]= 0; 1 <сШ.п; ++) (с ПИ <0)е = 
Лог (ши /= 0; 1 < с2.; ++) (с? [Л < 0) е2 [7] = 0;} 
ге ауегС1 = (пЕ2А(ТА)(с1) — с с1= с1 - ауегСТ; 
ге ауегС2 = (пе2А(Тй)(с2) — та5вас2)/агеа0; 2 = с2 — ауегС?2; 
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Визуализация расчетов. Визуализация расчетов на каждом шаге по времени обеспечи- 
валась кодом программы р/о!(яета, стт = +1) ‚ выводившим на экран линии уровня функции с, 


в полной мере отражающей распределение концентраций биополимеров. 

Особое внимание было уделено вопросу о задании начального распределения концентра- 
ций биополимеров, т.е. формы, размеров и расположения в ЭК пятна смеси разделяемых ве- 
ществ. Начальное распределение полимеров задавали либо в виде круга [5]: 


(у) А [1+8 (4 (<) +6-») -*)), 
(А, — амплитуда, В — параметр сглаживания, (Ур) — центр пятна, х, — его радиус), либо в ви- 
де прямоугольника со сглаженными углами: 
6 (х,у)= 0,254, | (В, (хх) (в, (х-х,)) [№ (в, (у- у, ))- (в, (-х,)) |, 
(В, В, - параметры сглаживания углов прямоугольника, занимающего область 
{(х,у):х <х ху <узу,}). 
Исследование динамики концентраций биополимеров. /7ример 1 (рис.3). Расчеты прово- 


дили для десяти гипотетических аминокислот с равномерным распределением изоэлектрических 
точек: Др/ = 0,3, ДрК = 0,15. При этом предполагалось, что иммобилизованный градиент задан 





линейной функцией. 





Т=0,1 


2 




















Рис.3. Процесс ИЭФ десяти аминокислот с равномерным распределением изоточек: 
Др! = 0,3 , АрК = 0,15 ; градиент РН линейный, рН, =7,0, рН, =10,0 (начало): 


а- первый этап; б- второй этап 
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Рис.3. Окончание 


Из рис.3 видно, как осуществляется процесс расслоения исходной смеси на фазы. Исход- 
ная круговая форма пятна со временем трансформируется в эллиптическую, а затем фракция 
приобретает вид вертикальной полосы. По окончании расслоения процесс приобретает стацио- 
нарный характер, так как для Т =3,5 и Т =4,0 формы пятен имеют минимальные различия. 

Численный эксперимент показал, что на итоговую картину не влияет исходная форма пят- 
на смеси — для распределения в виде прямоугольника наблюдалась картина, аналогичная 
рис.3. Результаты проведенного исследования подтверждаются результатами одномерного моде- 
лирования ИЭФ [6]. 

Пример 2 (рис.4). Проведены расчеты для восьми стандартных аминокислот А5р, 
м-АБК, а- Азр - Ня, Ня -Су-Г, НБ-СЬ-П,В-Аа- НБ, Ту"- Ате, От. 

В расчетах были использованы характеристики биополимеров-носителей [7] (см. таблицу). 
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Азр, м-АБК, 2-Йзр-Не, НИ-ВМ +, Н-СРу-Щ, 6 -АГа Ник, Туг-Аго, Огни: 
РН =2,5 рН>=105 





Т=0, 1 
























































Рис.4. Процесс ИЭФ восьми стандартных аминокислот А5р, м- АБК, &-— Азр- Ня, Нк -Су-Г, 
Но -Су-П ‚ В- Аа - Нв5 ‚ Туг - Ага , Отт ; градиент рН -— линейный, РН, =2,5, РН, =10,5 


Характеристики стандартных аминокислот 















































Амино- ве яко р Арк Коэффициент 

ВИСлОта подвижности, х10“ 
Ар 1,88 3,65 2,77 1,77 3,078 
м-АБК 3,12 4,74 3,93 1,62 3,119 
9— Азр- Н 3,02 6,82 4,92 3,80 2,187 
Н& - Су-1 5,28 7,8 6,80 2,02 2,487 
Ни -Су-П 6,27 8,57 7,42 2,30 2,487 
В- 4а- НБ 6,83 9,51 8,17 2,68 2,3837 
Туг — Ато 7,55 9,80 8,68 2,25 1,638 
Оти 8,65 10,76 9,70 2,11 3,223 
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Из рис.4 видно, что наиболее четкое разделение наблюдается для аминокислот А5бр, 
м-АБК и а- Азр- Н5, поскольку шаг по АрГ для них максимален. В то же время, 
В- Аа -Н15 и Туг- Аге выделены преимущественно в виде смеси из-за относительно близких 


изоэлектрических точек. Таким образом, созданное программное обеспечение позволяет предска- 
зать поведение сложных систем биополимеров в заданных условиях ИЭФ. 

Выводы. 1. Построена универсальная двумерная математическая модель ИЭФ в прямоугольной 
электрофоретической камере с фиксированным (иммобилизованным) градиентом рН, описываю- 
щая поведение в процессе ИЭФ заданного числа произвольных биополимеров, для которых долж- 
ны быть известны константы диссоциации рК“, рК‘®, а также коэффициенты миграции ц,. 


Кроме того, для работы с моделью необходимо задать градиент рН в ЭК. 

Математическая модель представляет собой краевую задачу, для решения которой были 
использованы средства интегрированной среды разработки (языка) ЕгееРет++. На основе 
стандартных средств языка созданы собственные алгоритмы решения и программа, визуализи- 
рующая результаты расчетов и позволяющая наблюдать за трансформацией пятен биополимеров 
во времени. 

Построенная модель позволяет исследовать динамику концентраций биополимеров 
в зависимости от параметров ЭК - разности потенциалов и градиента рН; параметров самих био- 
полимеров — констант диссоциации, коэффициентов миграции, а также начальной формы пятна 
смеси. 

2. Проведенные расчеты показали, что разделение аминокислот происходит тем 
лучше, чем меньше значения АрК’ по сравнению с длиной интервала между их изоэлектрически- 


ми точками. 

3. На основании расчетов можно прогнозировать результат ИЭФ для заданного режима ЭК 
и системы биополимеров. С помощью модели также можно оптимизировать эксперимент, варьи- 
руя состав электролита, разность потенциалов и рН в ЭК. Таким образом, построенная модель 
может быть использована для оптимизации реального эксперимента, повышения разрешающей 
способности метода, а также экономии времени и средств. 
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ТМ/О-ОТМЕМ$ТОМАЕ МАТНЕМАТТСАЕ МОБЕММС ОЕ ТЕР ВУ ЕгееРет+-+ 
1МТЕСКАТЕО ОЕУЕГОРМЕМТ ЕМУТВОММЕМТ 


Е.М. ЗАКНАКО\УА 
(Адтиа! УзНаКоу Майите Зае Асадету, Ко$ю\-оп-Ооп Бгапсй) 


Тре ргоМет оЁ 2)0-@тВивоп оЁ Бороутегз Бу пе ТЕР т те гечапди/аг еесоуйс се! ийй {пе тторИтеа 
рН-дгафепЕ 15 зо№еа. Тве сотриег ргодгат сгеае4 оп те Базе оЁ ЕгееРет++ регтйз ю ГпуезНдее {етрога/у 
{Пе 1ЕЕ-ргосез$ апа 0 тодЕ! апа орйтие а геа/ ехрейтепЕ 

Кеуигога5: ПЕР, та етайса/ зиптшайоп, Ппйе веетепЁ те од, тедгаеа деуеортепЁЕ епигоптепё, 
НгееРет++. 
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